package scu.maqiang.peridynamics.interpolation;

import scu.maqiang.mesh.Tecplot;
import scu.maqiang.numeric.Algorithms;
import scu.maqiang.numeric.DirectDSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.ParamCheck;

import java.util.Arrays;

public class PDApproximation {
    public static double[] polynomials(double[] x, int order) {
        double x1, x2, x3, x4;
        int dim = x.length;
        if(dim == 1) {
            switch (order) {
                case 0:
                    return new double[]{1.0};
                case 1:
                    x1 = x[0];
                    return new double[]{1.0, x1};
                case 2:
                    x1 = x[0];
                    return new double[]{1.0, x1, x1 * x1};
                case 3:
                    x1 = x[0];
                    x2 = x1 * x1;
                    x3 = x1 * x2;
                    return new double[]{1.0, x1, x2, x3};
                case 4:
                    x1 = x[0];
                    x2 = x1 * x1;
                    x3 = x1 * x2;
                    x4 = x2 * x2;
                    return new double[]{1.0, x1, x2, x3, x4};
                default:
                    throw new IllegalArgumentException("order should be from 0 to 4");
            }
        }
        if(dim == 2) {
            x1 = x[0];
            x2 = x[1];
            switch (order) {
                case 0:
                    return new double[]{1.0};
                case 1:
                    return new double[]{1.0, x1, x2};
                case 2:
                    return new double[]{1.0, x1, x2, x1 * x1, x1 * x2, x2 * x2};
                case 3:
                    return new double[]{1.0, x1, x2, x1 * x1, x1 * x2, x2 * x2,
                                       x1 * x1 * x1, x1 * x1 * x2, x1 * x2 * x2, x2 * x2 * x2};
                case 4:
                    return new double[]{1.0,
                                        x1, x2,
                                        x1 * x1, x1 * x2, x2 * x2,
                                        x1 * x1 * x1, x1 * x1 * x2, x1 * x2 * x2, x2 * x2 * x2,
                                        x1 * x1 * x1 * x1, x1 * x1 * x1 * x2, x1 * x1 * x2 * x2, x1 * x2 * x2 * x2, x2 * x2 * x2 * x2};
                default:
                    throw new IllegalArgumentException("order should be from 0 to 4");
            }
        }
        if (dim == 3) {
            x1 = x[0];
            x2 = x[1];
            x3 = x[2];
            switch (order) {
                case 0:
                    return new double[]{1.0};
                case 1:
                    return new double[]{1.0,
                                        x1, x2, x3};
                case 2:
                    return new double[]{1.0,
                                        x1, x2, x3,
                                        x1 * x1, x1 * x2, x1 * x3, x2 * x2, x2 * x3, x3 * x3};
                case 3:
                    return new double[]{1.0,
                                        x1, x2, x3,
                                        x1 * x1, x1 * x2, x1 * x3, x2 * x2, x2 * x3, x3 * x3,
                                        x1 * x1 * x1, x1 * x1 * x2, x1 * x1 * x3, x1 * x2 * x2, x1 * x2 * x3,
                                        x1 * x3 * x3, x2 * x2 * x2, x2 * x2 * x3, x2 * x3 * x3, x3 * x3 * x3};
                case 4:
                    return new double[]{1.0,
                                        x1, x2, x3,
                                        x1 * x1, x1 * x2, x1 * x3, x2 * x2, x2 * x3, x3 * x3,
                                        x1 * x1 * x1, x1 * x1 * x2, x1 * x1 * x3, x1 * x2 * x2, x1 * x2 * x3,
                                        x1 * x3 * x3, x2 * x2 * x2, x2 * x2 * x3, x2 * x3 * x3, x3 * x3 * x3,
                                        x1 * x1 * x1 * x1, x1 * x1 * x1 * x2, x1 * x1 * x1 * x3,
                                        x1 * x1 * x2 * x2, x1 * x1 * x2 * x3, x1 * x1 * x3 * x3,
                                        x1 * x2 * x2 * x2, x1 * x2 * x2 * x3, x1 * x2 * x3 * x3,
                                        x1 * x3 * x3 * x3, x2 * x2 * x2 * x2, x2 * x2 * x2 * x3,
                                        x2 * x2 * x3 * x3, x2 * x3 * x3 * x3, x3 * x3 * x3 * x3};
                default:
                    throw new IllegalArgumentException("order should be from 0 to 4");
            }
        }
        throw new IllegalArgumentException("Dimension should be from 1 to 3");
    }
    public static double[] Interpolation1D(double[] inputX, double[] inputY, double[] outX, int order) {
        int inputSize = inputX.length;
        //2. 输出参数
        int outN = outX.length;
        double outDx = (inputX[inputSize - 1] - inputX[0]) / outN;
        //3. 计算Lj_T
        double Eps = 1.0e-9;
        int p = 3;
        double[] Lj_T = new double[inputSize];
        double[] rho = new double[inputSize];
        for(int k = 0; k < outN; k++) {
            int kk = k;
            Arrays.setAll(rho, i -> 1.0 / (Math.pow(Math.abs(outX[kk] - inputX[i]), 3) + Eps));
            double rhoSum = MVO.sum(rho);
            Arrays.setAll(Lj_T, i -> Lj_T[i] + rho[i] / rhoSum * outDx);
        }

        double[] outY = new double[outN];
        double[] xi_kj = new double[inputSize];
        int[] distanceOrder = new int[inputSize];
        int numberOfFamily = inputSize;
        double delta_k = 0.0;
        int N = order + 1;
        double[][] A = new double[N][N];
        double[] b = new double[N];
        b[0] = 1.0;
        double[] P = new double[N];
        double[] weight = new double[N];
        double temp = 0.0;
        for(int k = 0; k < outN; k++) {
            int kk = k;
            MVO.fill(A, 0.0);
            Arrays.setAll(xi_kj, i->Math.abs(outX[kk] - inputX[i]));
            Algorithms.selectionSort(xi_kj, distanceOrder, true);
            delta_k = MVO.max(xi_kj);
            for(int jmem = 0; jmem < numberOfFamily; jmem++) {
                int j = distanceOrder[jmem];
                double xi = inputX[j] - outX[k];
                Arrays.setAll(P, i -> Math.pow(xi, i * 1.0));
                temp = delta_k /(Math.abs(xi) + 1.0e-6);
                Arrays.fill(weight, temp * temp);
                for(int ii = 0; ii < N; ii++) {
                    for(int jj = 0; jj < N; jj++) {
                        A[ii][jj] += P[ii] * P[jj] * weight[jj] * Lj_T[j];
                    }
                }
            }
            Arrays.fill(b, 0.0);
            b[0] = 1.0;
            DirectDSolver solver = new DirectDSolver(A);
            double[] Acoef = solver.PGaussSolve(b);
            for(int jmem = 0; jmem < numberOfFamily; jmem++) {
                int j = distanceOrder[jmem];
                double xi = inputX[j] - outX[k];
                double yVal = inputY[j];
                Arrays.setAll(P, i -> Math.pow(xi, i * 1.0));
                temp = delta_k /(Math.abs(xi) + 1.0e-6);
                Arrays.fill(weight, temp * temp);
                double[] tempVec = new double[P.length];
                Arrays.setAll(tempVec, i -> Acoef[i] * P[i] * weight[i]);
                outY[k] += yVal * MVO.sum(tempVec) * Lj_T[j];
            }
        }
        return outY;
    }

    public static double[] Interpolation2D(double[][] inputXY, double[] inputZ, double[][] outXY, int order) {
        return null;
    }

    public static double[] Interpolation3D(double[][] inputXY, double[] inputZ, double[][] outXY, int order) {
        return null;
    }

    public static double[] Fit1D(double[] inputX, double[] inputY, double[] outX, int order) {
        return null;
    }

    public static double[] Fit2D(double[][] inputXY, double[] inputZ, double[][] outXY, int order) {
        return null;
    }

    public static double[] Fit3D(double[][] inputXY, double[] inputZ, double[][] outXY, int order) {
        return null;
    }


    public static void main(String[] args) {
        //1. 输入参数
        double[] inputX = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
        double[] inputY = {2, 2, 2, 2, 2, 3.2, 3.2, 3.2, 3.2, 3.2};
        int inputSize = inputX.length;

        //2. 输出参数
        int outN = 200;
        double outDx = (inputX[inputSize - 1] - inputX[0]) / outN;
        double[] outX = MVO.linspace(0.5 * outDx, inputX[inputSize - 1] - 0.5 * outDx, outN);

        double[] outY2 = PDApproximation.Interpolation1D(inputX, inputY, outX, 2);
        double[] outY3 = PDApproximation.Interpolation1D(inputX, inputY, outX, 3);
        double[] outY4 = PDApproximation.Interpolation1D(inputX, inputY, outX, 4);
        double[] outY5 = PDApproximation.Interpolation1D(inputX, inputY, outX, 5);
        Tecplot.LineXY("NewInterpolation1D.dat", outX, new double[][]{outY2, outY3, outY4, outY5});
    }
}
